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ABSTRACT 

A formal derivation is presented of the energy transfer rate between radiation and 
matter due to the scattering of an isotropic distribution of resonant photons. The 
derivation is developed in the context of the two-level atom in the absence of collisions 
and radiative transitions to and from the continuum, but includes the full angle- 
averaged redistribution function for photon scattering. The result is compared with 
previous derivations, all of which have been based on the Fokker-Planck approxima- 
tion to the radiative transfer equation. A new Fokker-Planck approximation, including 
an extension to higher (post-diffusive) orders, is derived to solve the radiative transfer 
equation, and time-dependent numerical solutions are found. The relaxation of the 
colour temperature to the matter temperature is computed as the radiation field ap- 
proaches statistical equilibrium through scattering. The results are discussed in the 
context of the Wouthuysen-Field mechanism for coupling the 21cm spin temperature 
of neutral hydrogen to the kinetic temperature of the gas through Lya scattering. 
The evolution of the heating rate is also computed, and shown to diminish as the gas 
approaches statistical equilibrium. 

Key words: atomic processess - cosmology: theory - line: formation - radiative 
transfer - radio lines: general - scattering 



1 INTRODUCTION 



The prospects for detecting the Intergalactic Medium (IGM) 
prior to reionization by extremely large radio telescopes 
like the Low Frequency Array (LOFAR) and the Square 
Kilometre Array (SKA) has led to a renewed interest the 
Wouthuysen-Field mechanism for coupling the spin tem- 
perature of neutral hydrogen to the kinetic temperature of 
the gas through the scattering of Lya photons (Wouthuy- 
sen 1952; Field 1958). Following Field's (Field 1959a, b), 
analysis of the scattering of Lya photons in a homogeneous 
medium and their role in setting the spin temperature of 
neutral hydrogen (Field 1958), the physics of 21cm exci- 
tation and absorption was re-examined in the light of cur- 
rent cosmological models by Madau, Meiksin, & Rees (1997, 
hereafter MMR) and Tozzi et al. (2000). A large number 
of various scenarios and effects have since been explored in 
the literature (see Bowman, Morales & Hewett 2005 and 
references therein). 

In addition to de-coupling the hydrogen spin tempera- 
ture from the temperature of the Cosmic Microwave Back- 
ground (CMB), MMR showed that Lya photons were in 
principle capable of heating cold hydrogen gas through re- 
coils. The derivation was based on the momentum trans- 
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ferred to the hydrogen atoms by the Lya photons after re- 
coil (Field 1959b, Basko 1981). MMR supported their ar- 
gument by identifying the recoil heating term in the ra- 
diative transfer equation for the scattering of the photons 
using a Fokker-Planck approximation with atomic recoil 
of Rybicki & Dell' Antonio (1994), who extended previous 
Fokker-Planck treatments (Unno 1955; Harrington 1973; 
Basko 1981; Krolik 1990) to include the full Voigt line pro- 
file. While recognising that the energy exchange would cease 
once the radiation and matter reached thermal equilibrium, 
MMR did not discuss the role of photon diffusion through 
frequency on the energy exchange rate, as its effect on energy 
exchange is negligible for sufficiently cold gas. As the radi- 
ation field approaches statistical equilibrium through scat- 
terings, however, the photon diffusion term reduces the net 
heating rate, which ultimately vanishes once the radiation 
field achieves statistical equilibrium at frequencies near res- 
onance, with a colour temperature that matches the matter 
temperature. This will generally occur well before the radia- 
tion field reaches thermal equilibrium with the matter. The 
degree to which the approach to statistical equilibrium re- 
duces the heating rate depends on the frequency range over 
which equilibrium is achieved. While it may be achieved in 
the core of the absorption profile on a timescale of the or- 
der of the mean free time for scattering at line-centre, it will 
take longer outside the core and significant heating may still 
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persist. A central purpose of this paper is to derive the en- 
ergy exchange rate from the full radiative transfer equation, 
justify the use of the Fokker-Planck equation for describing 
the exchange of energy between the radiation and the mat- 
ter, and compute the evolution of the heating rate in this 
approximation for some simple applications. This is done in 
the approximation of the two-level atom, neglecting colli- 
sions and scatterings to and from the continuum. A second 
purpose it to compute the timescale over which the colour 
temperature relaxes to the kinetic temperature of the gas, as 
this is crucial for determining the signature of 21cm emission 
or absorption through the Wouthuysen-Field mechanism. 

Under the assumption that the Fokker-Planck equation 
of Rybicki & Dell' Antonio (1994) adequately described the 
exchange of energy with the matter, Chen & Miralda-Escude 
(2004) applied Rybicki & Dell'Antonio's cosmological ver- 
sion of the Fokker-Planck equation to the heating of the 
IGM, assuming that the comoving radiation density had 
reached a steady-state. In this case, the energy exchange 
between the radiation field and the matter virtually ceases, 
as in the static case, except for a residual amount arising as 
photons are redshifted through the resonance frequency. No 
estimate, however, is made of the time required to reach the 
equilibrium heating rate. 

Following Deguchi & Watson's (1985) formulation of 
the resonant photon scattering problem including stimulated 
emission, Rybicki (2006) presents an alternative derivation 
of the energy exchange rate using a Fokker-Planck approx- 
imation, and derives an approximate form for the radiative 
transfer equation that is structurally similar to the Kompa- 
neets equation. 

In order to solve the radiative transfer equation, a dif- 
ferent Fokker-Planck approach is adopted here. The Fokker- 
Planck approximation derived is more complete than any 
previously formulated in this context. In particular, the for- 
mulation presented here explicitly conserves photon num- 
ber and energy to second order in derivatives of the radia- 
tion field, as well as when extended to higher orders, per- 
mitting more accurate solutions to be obtained in a self- 
consistent manner than was previously possible. This reme- 
dies a shortcoming of the Fokker-Planck treatment of Ry- 
bicki & DelP Antonio (1994), for which particle conserva- 
tion was imposed by hand. Since then Rybicki (2006) t has 
presented an improved Fokker-Planck treatment. There are, 
however, some differences between the Fokker-Planck treat- 
ment provided here and Rybicki's. These are discussed be- 
low. 

In the next section, the basic theory of resonance scat- 
tering is reviewed, including a derivation of the rate of en- 
ergy exchange between the radiation field and the scattering 
medium. The Fokker-Planck approximation is discussed in 
§3, and applied to some simple problems in §4. A discussion 
of the results and conclusions are presented in §5. 



2 BASIC THEORY 

2.1 Radiative transfer and heating equations 

Statements of the radiative transfer equation in the litera- 
ture are generally approximations, and not all are adequate 
for formulating the full energy transfer problem between the 
photons and the scattering medium. A standard form of the 
radiative transfer equation in terms of the photon energy 
density (Mihalas 1978, §2) in fact leads to the complete ab- 
sence of energy transfer due to the implicit approximations 
made. In this section, a derivation of the radiative trans- 
fer equation satisfactory for treating the effects of energy 
exchange by the scatter of resonant photons is presented. 
A formal expression for the energy exchange rate for an 
isotropic radiation field is derived in this section which uses 
the full scattering integral, unlike all previous discussions 
which were based on the Fokker-Planck approximation. 

In order to clarify the assumptions made in the radiative 
transfer equation for resonant photons, a derivation from 
basic principles is provided. The problem is formulated in a 
non-relativistic context. Polarisation is suppressed as well, 
but the results are readily generalised. 

In a static medium, the rate of change of the energy 
in the radiation field at frequency v at position r at time t 
moving (in the laboratory frame) into direction fi per unit 
area per unit solid angle per unit time per unit frequency is 
(Mihalas 1978, Section 2.2) 



Dt 



dly(r,t,n) 
dt 



+ cVIu(r,t,n) 



= C7/„(r, n,t) -cx v (T,n,t)I v (r,t, n), (1) 

where I v (r, t, n) is the specific intensity, rj v is the specific 
emissivity (energy emitted per unit volume per unit solid 
angle per unit time per unit frequency), Xv is the spe- 
cific extinction (fraction of incident energy removed per unit 
length), and c is the speed of light. (Finite light travel-time 
effects are neglected.) 

For simplicity, it will be assumed that the emissivity and 
extinction both arise entirely due to scattering by a two-level 
atom. (Much of the following discussion is based on Mihalas 
1978, §13-4. See also Milkey & Mihalas 1973.) Scatterings 
to and from the continuum levels are neglected, as are ex- 
citations and ionizations due to collisions. Local isotropy of 
the radiation field is also assumed, although global varia- 
tions with position are allowed. Under these assumptions, 
the rate of change of atoms in the upper state due to inter- 
actions involving the production of photons of energy hu is 
given by 

dn u (v) 



dt 



= -n u (z/) (A u i + B u i Ju) 

p oo 

+ niBi u / dv R(y' ,v)J v > , 
Jo 



(2) 



t The present paper was substantially complete, but not submit- 
ted, prior to the posting of Rybicki (2006). Differences between 
the present work and Rybicki's are commented on in the text. 



where J„ — <f> (duj /4tt) I„ is the usual angle-averaged inten- 
sity. (J„ = I v for an isotropic radiation field.) The terms 
on the right hand side correspond, respectively, to sponta- 
neous decays given by the Einstein coefficient A u i, stimu- 
lated emissions described by the downward coefficient B u i 
and scatterings from the lower to upper level described by 
the upward coefficient Bi u . The scattering function R(y' , v) 
is the angle- averaged probability that a photon of frequency 
v' is scattered to frequency v by an atom in the lower state. 
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It is normalised by ip(v) = J"°° dv ' R(v, v'), where ip(v) is 
the absorption profile of the atoms in the lower state. The 
function n u {v)dv describes the number density of atoms 
in the upper state capable of emitting photons in the fre- 
quency range (v, v + dv) in the observer's frame. In general 
n u (v)/n u , where n u is the total number density of photons 
in the upper level, is not given by the natural emission profile 
if>(v) — f Q dv' R(v' ,v) and must be solved for. 
The emissivity is then given by 



M 4^' 



r\ v = hvn 
and the absorption coefficient by 
hv 



(3) 



X, 



Biu . . B u i 



(4) 



Expressed in terms of the radiation energy density u v 
4nrJ v /c, the radiative transfer equation becomes 



Du v (v, t) 
Dt 



duvir, t) 



+ cV«„(r, t) 



dt 

hvn u (v)A u i 

B ut n u (v) 



hvni—^- ctp(v)u u 
4ir 



1 - 



Biu ni(p(v) 



(5) 



In most applications, the level populations will have 
reached equilibrium. Eq. (2) then gives 



Bin f dv' R(v', v)u v i 



n u (v)=m , 

4ttA u i/c + B ui u v 

Noting that the photon occupation number is M(y 

e 'u v i '(8nhv 3 ] 
as 

Du v (v, t) 



(6) 



, the rate equation for u v may be expressed 
du„ (r, t) 



Dt 



dt 

hvni(v, t)c 



+ cVw„(r, t) 
B lu f dv'R(v',v)u„,(r,t) 



(7) 



4tt l+M{v,r,t) 



~hvni (r, t)c?^-ip(v)u v 
/ dv' R(v',v) 



N{v',v,t) 



<p(v) [1+M(v,r,t)} 



after using the relations between the Einstein coefficients. 
When the radiation reaches thermodynamic equilibrium, it 
becomes blackbody radiation at the same temperature as 
the matter, so that Du„/Dt = 0. This permits a relation to 
be derived between R(v, v') and R(v' , v) in detailed balance. 
Setting Af(v) = l/[exp(hv/kT) - 1] gives 

3 hv' /kT 



R(v',v) 
R{v, v') 

where the limit applies for hv » kT and hv' » kT cor- 
responding to the classical limit of a dilute radiation field. 
This limiting form corresponds to the result of Deguchi & 
Watson (1985) (noting the change in their normalization 
of R(v,v')). They claim their result is general, although it 
seems to apply only for low values of the photon occupancy. 

Eq. (7) simplifies. Defining the coefficient Xo(r, t) = 
ni(r,t)(Bi u /4:iT)hvo, where vq is the resonant frequency, af- 



ter simplification the equation of radiative transfer becomes 



Du v {v,t) du„{r,t) 
~~ dt 
Xo(r,i)c 



Dt 



hv 
-ip(v)u v (r,t) 



+ cV-u„(r, t) 
hv 



dv R{y' , v)u v i 
+ hvS{v). (9) 



A source term S(v) of photons (per unit volume) has also 
been added (which may be taken to be an arbitrary func- 
tion over space and time). It is noted that in the presence 
of additional interactions that are capable of redistributing 
the level populations, such as collisions or ionizations and 
recombinations to and from the continuum, the relations 
become considerably more involved. In general it is neces- 
sary to solve for the detailed level populations including the 
emission profile in order to establish the fraction of atoms 
in the lower and upper states. (See Mihalas 1978 §13-4 for 
details.) 

An alternative form of the radiative transfer equation 
in terms of the number density of photons n„ = u v /hv is 
given by dividing Eq. (9) by hv 



Dn v (r, t) 
Dt 



dn v (v, t) 

dt 
Xo(r,t)c 
hv 



+ cVn„(r, t) 

/ dv' R(y' ,v)u v 
Jo 



-<fi(v)u„(r,t) 



+ S(v). 



(10) 



The rate at which energy is transferred to the matter 
per unit volume is given by integrating Eq. (9) over fre- 
quency 

D_ 

~bl 



G(r,t) = 



dvhvS(v) 



dvu v (v, t) 



Xo(r,t)c 



hv 



dv 



dv'hvR(v' , v)u v i (r, t) 



hv 



dvhvLp(y)u v (r, t) 



POO f-O 

/ dv \ 

JO JO 



dv (hv — hv)R(v' , v)u v f( 



after noting that J dvR(y' , v) — ip(v'). The final form shows 
that the energy transfer rate to the scattering medium is 
given by the rate at which the photon energies in the ob- 
server's frame shift during scattering. 

Eqs. (9) and (11) have some noteworthy symmetry 
properties that may be used to gain insight into the origin 
of the energy exchange in special situations. For resonant 
scattering, the redistribution function R(v',v) for most ap- 
plications is well-approximated as a function of (y — v'). 
Similarly tp(v) is an (even) function of (v — vo). For a 
source function S v that is a function of (v — vo), in the 
absence of an initial radiation field, the form of Eq. (9) 
shows that the radiation energy density u v may also be ex- 
pressed as a function u(v — vo). It is then possible to ex- 
press R(v, v') and u(y — vo) as sums of their even and odd 
parts in (y - v ), R(v,v') = R (e) (v,v') + R (o \v,v') and 
u(v-v ) =u {e) (v-vo) +u < - o) (v - v ), where R (e) (v,v') = 
[R{v,v') + R{v',v)]/2, R(°\v,v') = [R(v,v') - R{v',v)]/2, 
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tt (e '(V — vq) = [u{v — vq) + u(yo — v)]/2 and u^iy — vq) = 
[u(v— vq)— u{vq — f)]/2. The form of Eq. (11) then shows that 
energy exchange between the radiation field and the matter 
is driven only by the even-odd combinations R^u^ and 
/j>( o ) M 0) j n particular, if the source function Sv is an even 
function of (y — vq), then u(y — vo) will have an odd part 
only if R{u,v') does. If instead R(u, v') = R(v' ,v), as holds 
for the redistribution Cases I, II and III (Mihalas 1978, §13- 
3) without atomic recoil t, then u(y — vo) will be even in 
(y — vq) and no energy exchange will occur. It is evident from 
Eq. (8) that the asymmetry induced by atomic recoil gives 
rise to the energy exchange in this case. These symmetry 
relations are exploited below. 

It is shown in the Appendix that under the assumption 
of Case Ilr redistribution (perfectly sharp lower state, finite 
radiative decay time in upper state, frequency coherence in 
the restframe of the atom) , which is the case of most interest 
here, the integral over the redistribution function may be 
expressed in terms of the Voigt line profile ipv(u), so that 



G (r,t) = *°£l^ft 



f 

Jo 



dvu v {v, t) 



1 2 dipvM 

LJ — - 

2 dv 



(12) 



where the effects of atomic recoil have been included to order 
e = hvo/(2kTm a c 2 ) 1 '" 2 , where T is the gas temperature and 
m a is the mass of the recoiling atom, and co = vob/c, where 
b = (2kT/m a ) 1 ^ 2 is the Doppler parameter. Because tpv(v) 
is an even function of [y — v a ), Eq. (12) shows that energy 
transfer by recoil (the first term in brackets) is driven by the 
even part of u v = u(y — vo), while energy transfer through 
the second term is driven by the odd part. The odd part 
may also contain a term of order e, as it will for an even 
source. 

It is shown below that in the Fokkcr-Planck approxi- 
mation, a source function S v even in (y — v ) (or one well- 
approximated as even, such as a function very smoothly 
varying across the resonance frequency), will generate an 
odd part to u(y — vq) only through recoils. Thus both terms 
in Eq. (12) will be proportional to the recoil parameter e. 
It is also noted that any non-negligible variation in fre- 
quency of the continuum source across the resonant fre- 
quency may induce non-negligible energy transfer between 
the radiation and matter even in the absence of recoils. Ex- 
panding S v ?s S VQ + [y — h>o)(dS/dv)o shows that the correc- 
tion term [y — t>o)(dS/dv)o will drive an odd contribution to 
u v — u(y — vi$) in Eq. (9) even for a symmetric redistribution 
function (R[v, v'\ = R\y' ' , v\), and so create a non-vanishing 
energy exchange rate as described by Eq. (11). 

Following Rybicki (2006), it is useful to recast Eq. (12) 
in terms of the radiation colour temperature. Motivated by 
Field's (1958) introduction of an effective colour tempera- 
ture of the radiation field, MMR presented a formal defini- 
tion of colour temperature TV (called T a in their paper) in 
terms of the photon occupation number J\f{v): 



t Technically, redistribution Cases I, II and III do not allow for 
atomic recoils. The modifications with recoils shall be designated 
by adding the suffix V, thus, Case Ilr. 



T\r = —- 



d\ogN{v) 



dv 



(13) 



In order to account for stimulated scattering, Rybicki (2006) 
modified the temperature definition to T]^(y) = Tj^(v)[l + 
Miy)]. 

In light of the form of Eq. (12), it is now suggested 
that, in the absence of stimulated emission, the most natural 
definition for the colour temperature for describing energy 
exchange is in terms of the energy density of the radiation 
field. Through straightforward manipulations, Eq. (12) may 
be recast in the form 



G = P l n l ^hv I 1- 



where 



d Bin 



dvipv(y)u v 



(14) 



(15) 



is the total photon scattering rate per atom in the lower 
state, and the mean harmonic colour temperature has been 
defined by: 



(T U ) H = 



where 



T u (u) 



poo if °° 1 

duu v ip v {v) I duu v ip v {u) T — , 



d log u v 
dv 



(16) 



(17) 



The limit T « (T u )h of Eq. (14) recovers the heating rate 
given by MMR. At lower colour temperatures, the heating 
rate is suppressed by the efficiency factor 1 — T / (T U ) H . 

The temperature in Eq. (16) is identical to the "light 
temperature" Tl defined by Field (1958) in the context of 
the Wouthuysen-Field mechanism, as will now be demon- 
strated. If the total transition rates induced by the scat- 
tering of Lya photons from the lower hydrogen hyperfine 
n = 1 level (0S1/2) in hydrogen to the upper hyperfine level 
(15*1/2) is designated P^, and the reverse rate designated 
P£ , then Field (1958) (see also Tozzi et al. 2000), shows 
that in equilibrium, 



pL " ' " 



(«„!> + («„'}' 



(18) 



where {u Vi ) = J dvipv(y — Vi)u v is the energy density of 
the radiation field averaged over a Voigt profile centred 
on v = Vi. The four frequencies uo, v' , v\ and v[ cor- 
respond to the Lya transitions between the n = 1 and 
n = 2 pairs (o<Si/ 2 , 1 P1/2), (oS 1/2 ,i P3/2), (iS'i/2,1 P1/2), and 
(i»S'i/2)i T3/2)) respectively. Field (1958) defines the effective 
colour temperature Tl of the radiation field by P^ / P^o = 
3(1— T*/Tl), where fcT* = hvoi and hvoi is the energy differ- 
ence between the n = 1 hyperfine levels. A straightforward 
expansion of the Voigt profiles about ipv(v — fo) to lowest 
order in dipv/dv shows that Tl is identical to {T u )h, fur- 
ther justifying this as the most pragmatic definition of the 
colour temperature in the absence of stimulated emission (cf 
Deguchi & Watson 1985). 

Rybicki (2006) expresses the heating rate in a form sim- 
ilar to Eq. (14), except in terms of Tp(u) (his eq. [32]), in- 
stead of T u {v). He then approximates the scattering cross- 
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section as a 5-function to evaluate the colour temperature 
at line-centre to obtain (his eq. [35]): 



GRyb 



hv , 

ni ^hvo 

m a c 2 



(19) 



(In Rybicki's formulation, u v is approximated as hvon v in 
the frequency-integral in the definition of Pi, and stimu- 
lated emission is accounted for by multiplying the result by 
[1 + Af(v)], the combined effect of which is defined here as 
Pi* .) The derivation here shows that the result is more gen- 
eral, and in particular does not rely on the Fokker-Planck 
formulation. It does, however, depend on the specific form 
Eq. (12). 

The various definitions of colour temperature have rel- 
ative differences from T u of kT u /hv. Provided kT u « hv, 
the differences in these definitions are therefore small. More- 
over, because multiplicative factors of v were dropped in 
adopting the Voigt profile (so, for example, the Rayleigh 
scattering limit is not recovered), the actual multiplicative 
factors of v in defining the radiation colour temperature are 
not accounted for in this formalism. For these reasons, the 
approximation u v — hvon v is adopted below for the defini- 
tion of the colour temperature for convenience, to give T n (v) 
in analogy to Eq. (17): 



(Tn) H 

where 



f 

Jo 




dvn v ipv{v) I I dvn v tpv(v)- 

/() inxy) 



1 lm : 



(20) 



(21) 



2.2 Uniform expansion 

In a homogeneous and isotropic expanding medium, Eq. (9) 
must be modified by adding expansion and redshifting terms 
to duv/dt. For this case, Duv/Dt becomes (Peebles 1968; 
Peebles 1993), 



du v R 

dt ~ dt +6 r u,/ 



R du v 

'R~bV 



(22) 



where R(t) is the expansion factor. The equation may be 
simplified further by defining u v = (R/ Rq) 3 u v and S„ = 
(R/Ro) S v , where Ro is the expansion factor at some time 
to. In this case, the form of Eq. (9) is recovered if u v and S v 
are everywhere replaced by u v and S v , and Du v /Dt is now 
identified with 



Du v 
~Dt 



du v 



R du v 
"~R~l)v~ 



(23) 



A similar result follows for the number density n v by first 
expressing u v — hvn v in Eq. (10), and then replacing n„ 
everwhere by — (R/Ro) by Sv = (R/R ) 2 S V , and 

identifying Dnv/Dt with 



Dn v 
~Dt 



dh v 

~~dT 



R dh v 
V ~R~dv' 



(24) 



An integration of the rate equation for u v exactly recov- 
ers the heating rate of Eq. (11), after expressing Du v / Dt by 
Eq. (22), which allows for the energy lost by the radiation 
due to the expansion. 



2.3 Approximate formulations 

Several approximate form of the transfer equation for reso- 
nant radiation exist in the literature. A common form starts 
from Eq. (9), but sets hv = hvo in the coefficient, ignoring 
the frequency shift of a photon on scattering (Mihalas 1978, 
§2-1). The equation for the energy density u v then becomes 



Duy(r,t) 
Dt 



Xo(r,*)c 



f 

Jo 



dv ' R(v' , v)u v i (r, t) — (p(h>)uv{v, t) 



While adequate for describing the evolution of u v for most 
cases, this approximate form is inadequate for the purpose 
of computing the rate of energy transfer between the pho- 
tons and the scattering medium, since the rate should be 
proportional to the frequency shift per scattering event. 
A straightforward integration over frequency, noting that 
f dvR(v' ' ,v) — ip{v'), shows that this form for the scat- 
tering equation gives a vanishing value for the net energy 
exchange rate. 

An intermediate approximation between the exact form 
and the approximation above is provided by Basko (1978) 
for the photon number density n v . This starts from Eq. (10), 
but makes the approximations u v > = hvoTiv' and u v = hvon v 
on the right hand side to arrive at 



Dnv(r, t) 
Dt 



= Xo(r,t)c 



f 

Jo 



dv 1 R{v , v)n v i (r, t) — ip(i>)nv(r, t) 



Most Fokker-Planck treatments in the literature start with 
this form. It has the advantage over Eq. (25) of recovering an 
approximate form for the heating rate, given by multiplying 
by hv and integrating over frequency 



G(r,t) 



" f 

Jo 



dvS(v) 



D_ 

Dt 



dviiv(r, t) 



(2' 



= Xo(r,t)c I dv 
I a 



/ dv' (hv' 
Jo 



hv)R(v' , v)n v '(r, 



It is identical to Eq. (11) except that u v > has been approx- 
imated by hvon v > . Since the expression already allows for 
a frequency shift of the scattered photon, however, this ap- 
proximation introduces only a small, and normally negligi- 
ble, error. 



3 FOKKER-PLANCK APPROXIMATION 
3.1 Basic equation 

The integro-differential character of the radiative transfer 
equation with the full redistribution function makes its so- 
lution nearly intractable, requiring special numerical meth- 
ods. A useful technique is to obtain solutions in a Fokker- 
Planck approximation. Basko (1978) introduced a Fokker- 
Planck approximation for Eq. (26) for an infinite homo- 
geneous medium, although he included only the dominant 
contribution to the line profile in the wings. Krolik (1990) 
presents a formal expression for the Fokker-Planck approx- 
imation, but similarly evaluates the coefficients only in the 
wings. Rybicki & DelP Antonio generalised Basko's form by 
rederiving a Fokker-Planck equation retaining the full Voigt 
line profile. Following Unno (1955) and Harrington (1973), 
they expand n v > in Eq. (26) in a Taylor series about v = v' . 
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Their derived coefficients, however, did not explicitly con- 
serve photon number, and so they were forced to adjust the 
coefficients to ensure the number of photons is conserved. 
An alternative derivation was presented by Rybicki (2006) 
invoking detailed balance and stimulated emission. It is still, 
however, based on a Taylor series expansion of n v . 

A derivation of the Fokker-Planck equation is presented 
here which explicitly conserves both particle number and en- 
ergy while still retaining the full line profile. For the present 
formulation, it is useful to introduce the probability W(y, Q) 
for scattering a photon of frequency v to one of frequency 
v' = v — Q. Eq. (9) may then be re-expressed as 



Du v (r, i) 
Dt 



Xo(r, t) c 
hv 



hv 



}'W{v\Q')u v , 



f 

J — < 



dQW{v,Q)u v 



+ hvS(v), 



(28) 



where Q' = v' — v, and ip(v) — J dQW(v,Q) has been used. 
It is also assumed here and elsewhere that W (v, Q) becomes 
vanishingly small for large Q, so that the lower integration 
bound v' = in the first integral may be well-approximated 
by Q' —* — co. Unlike previous approximations, the standard 
Fokker-Planck approximation is based on a Taylor series ex- 
pansion of the product of both the scattering probability 
function and the distribution function. A Taylor series ex- 
pansion of W(v' , Q')u v i about v' = v gives, to second order 
in v' — v. 

v d 



W{v',Q')u v , w W{v,Q')u„ + Q'f- [W(v,Q')u„ 



dv 



(29) 



Substituting this in Eq. (28) gives the canonical Fokker- 
Planck form 



Duv(y, t) 
Dt 

where 



^^hv^{ {Q)r{")<<, (y.t) C50) 



+ 



~ [(Q 2 Mv)u v {v,t)]^+hvS{v), 



f 

J — f 



dQQ n W(v,Q) 



(31) 



and the dummy integration variable Q' has been replaced 
by Q. This form conserves energy to the accuracy of the 
approximation and to order b/c, as may be shown by an 
integration by parts over v and noting that W(v, Q) vanishes 
at the integration boundaries, 



Xo(r,t)c 



(32) 



Xo(r 



dvu u (v, t) 

dvh(Q)Lp(u)u v (v, t) 

t)c f°° f 00 

I dv I dv ' (hv — hv')R(y,v')u v (v,€) 

'o Jo Jo 



hv 

identical to Eq. (11). This justifies estimating the energy 
exchange rate between the photons and the gas using the 
Fokker-Planck approximation, provided the approximation 
yields an accurate solution for the radiation field. 



Dividing Eq. (30) through by hv gives the Fokker- 
Planck equation for the photon number density 



Dn„(v, t) 
Dt 



Xo(r,i)c d 



hv 
1 



dv 



(Q)<p{v)u„ 



(33) 



2 dv 



[<Q 2 M^K] } 



+ S(v), 



which is in explicit photon number conservation form (cf 
Krolik 1990). The coefficients (Q) and (Q 2 ) are evaluated in 
the Appendix for Case Ilr redistribution. For e = 0, they are 
identical to those derived by Rybicki & Dell' Antonio (1994) 
(noting that they expand in — Q in the notation here). They 
did not explicitly include the recoil terms, instead adding 
by hand the result from Basko (1981). Here the recoil terms 
are included self-consistently, resulting in additional higher 
order corrections. 

It is convenient to solve the equation in non-dimensional 
form. Denning the dimensionless frequency shift x = (v — 
vo)/co and the conformal time r = J Q dtxo(r, t)c/uj, and sub- 
stituting in the moments of Q in terms of the dimensionless 
Voigt profile 4>v{x) = (pv(y)w, Eq. (33) becomes 



Dn x (r,T) 
Dt 



D log LJ 
Dt 



n x (r, r) 



(34) 



d_ 

dx 



t(j> v {x) - 



t>v(x 



dx 



n x (r,t) 



1 8 2 
+ 2dx 2 



4 dcj>v(x) 

* v{x) -r^x- 



+ 

2 dx 



1 dtfyjx) 

3 dx 2 



n x (r,t) } +S(x) 
c) 



2e<j> v (x)n x (r,t) + 4> v (x)—n x (r,t) 



6dx 2 \ 



^ d<t> v (x) d cj> v (x) 



dx 



dx 2 



n x (r,t) 



+S(x) 

ld_ 

2 dx 



2e(/> v (x)n x (r, t) + <j> v (x)—n x (v, t) 



+ 



2 , 

+ r 

6 dx 
-4e 



2 <j>v(x) 
dx 2 



n x (r,t) + 



ld 2 4> v (x) d 



6c d 2 <t> v (x) 
dx 2 



3 dx 2 

d 3 (p v (x) 
dx 3 



n x (r, t) 



dx 
n x (r,t) 



fivjx) d 
dx dx 



n x (r,t) }+S(x), 



where ti x (v,t) = n v {v, t)u> and S(x) = S(v)u) /(xoc) have 
been denned, a possible time-evolution of co has been allowed 
for, and terms of order xb/c have been neglected, consistent 
with the adoption of the Voigt profile (see the Appendix). 
It should be noted that the valid recovery of the effects of 
atomic recoil in this approximation requires b/c « 2e. For 
hydrogen Lya , this corresponds to T << hvo/k ~ 10 K. 

The approximations made by Rybicki & DellAntonio 
(1994) correspond to excluding the term in curly braces in 
the second equality of Eq. (34). The result has the form 
of a diffusion equation with diffusion coefficient <f>v(x)/2. 
The same form results from the full Fokker-Planck equation 
when higher order derivatives of <f>v{x) are neglected, as is 
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appropriate far in the wings of the profile (Harrington f973; 
Basko f978, 1981). It will be referred to here as the Fokker- 
Planck diffusion approximation (FPDA), to distinguish it 
from the standard Fokker-Planck approximation (FPA). 

The formulation of Rybicki (2006) (his eq. [20]), corre- 
sponds to excluding the term in curly braces in the third 
equality of Eq. (34). It is another FPDA equation. It is 
adopted here with the diffusion coefficient D(x) = <pv(x) + 
(l/3)d 2 cf>v(x)/dx 2 (his eq. [11]). The equation also includes 
an extra term of order b/c (arising from the term —2n v jv in 
his eq. [20]), and a quadratic term in n x arising from stimu- 
lated emission. These latter two terms, not included in the 
discussion here, play pivotal roles when the gas tempera- 
ture approaches the equivalent temperature of the resonant 
frequency or as the radiation thermalises with the matter, 
which, in the case of an optically thick medium, results in 
blackbody radiation. Comparisons with the approximation 
of Rybicki (2006) are confined here to the diffusion equation 
terms using the diffusion coefficient above. 

3.2 Higher order extension 

The above formalism is readily extended to higher orders. 
This is useful when solving the equations to test the degree 
of convergence of the original second-order Fokker-Planck 
approximation. Extended to fourth order (by continuing the 
expansion in Eq. (29) to fourth order), Eq. (34) takes the 
form 

Dn x (v,r) Dlogoj 



Dt 



Dt 



-n x (r,r) 



(35) 



d_ 

dx 



1 dcj>v(x) 
^ x) -2-^x~ 



n x (r, t) 



+ - 



2dx 2 



, . 4 d(j) V {x) 



3 dx 2 



n x (r,t) 



+ - 



6 da; 3 
1 d<&(x) 



4 dx 3 



e<f)v{x) - 
n x (r,t) 



^ d(f> v {x) 3^ d<j)y(x) 
dx 2 dx 2 



+ 7T7 



1 

24 dx 4 
8 d(j>lr(x) 



#v(a;) _l2e^M + 3- d ^ (a;) 



5 e " 



dx s 



5 dx 4 



dx 
n x (r,t) 



dx 2 



+ S(x) 



The third and fourth order terms correspond, respectively, 
to the terms with third and fourth derivates of the expres- 
sions in curly braces. Explicit expressions are provided in 
the Appendix for the moments {Q 3 } and (Q 4 ) from which 
the higher order terms above are derived. 

3.3 Uniform expansion 

Eqs. (34) and (35) continue to apply in an isotropic and 
homogeneous expanding medium if n x is replaced by h x = 
(R/Ro) n x and Dn x /Dr is given by 



Dn x (r) dn x (r) dn x (r) 
- 7" 



Dr 
where 



dr 



dx 



(36) 



7 = 



v H 8tt H 



Xoc 3 \lniA ul ' 



(37) 



is the ratio of the scattering time at line-centre to the ex- 
pansion time H and is known as the Sobolev parameter 
(Chugai 1980; Rybicki & Dell' Antonio 1994). It is also the 
inverse of the optical depth at line-centre Ao in a homoge- 
neous and isotropic expanding universe with Hubble param- 
eter H = R/R (Field 1959a). 



3.4 Method of solution 

It is instructive to express Eq. (35) in Fourier space. Denot- 
ing the Fourier transform of n x (r, r) by h K (r, r), Eq. (35) 
becomes 



Dh K (r, t) 
Dt 



D log uj „ . . 

-n K (r, r) = 



Dt 



(38) 



K 

2^ 



da 



+ i r 2 a(a 2 -8) + l K 3 (4-3a 2 + ia 4 ) 

1 2 / 2 8\ 

r r -3) 



l + -na+ 4 . 



1 3 / 2 2 \ 



'{a)h K - a {r,T), 



where 4>v(a) is the Fourier Transform of 4>v{x). Each or- 
der of n corresponds to the same order in the Taylor series 
expansion. The difference between the second order form 
above for the non-recoil terms and the Rybicki (2006) form 
is that in the latter a/2 above is replaced by a/2 — a A /6. 

The equation is solved by fourth-order Runge-Kutta 
with an adaptive timestep. The scheme is not fully stable for 
terms beyond the third order. While adequate for producing 
the solutions presented in this paper, it is not sufficiently ef- 
ficient to integrate to much longer times. Ultimately a mod- 
ified Crank-Nicholson method would be most desirable for 
solving Eq. (35). Such an approach is under investigation. 

The heating rate expressed in Fourier space may be de- 
rived from Eq. (12) and noting that 4>v{ — a) = 4>v{k), 

G(r,t) = f dn (e-ii«)0v(K)Mr,r) (39) 



2tt 

Xoch 
2ty 



dK 



eni r) (r, r) + ^/tn^(r,r) 



- (0 1 



where ni r ' and fi^ are the real and imaginary parts of h K , 
noting that the real and imaginary parts are, respectively, 
even and odd functions of k since n x is real. The heating 
rate may also be evaluated using Eq. (14). To the accuracy of 
approximation in the treatment here, (T u )h may be replaced 
by (Tn)H- Expressed as a ratio of integrals in Fourier space, 



{T n ) H = - 



k dK ' k4>v (k) 



(40) 



Cosmological expansion is accounted for by identifying 
Dh K /DT with 



DnJj) _ dnjj) 
Dt ~ Ot 



+ i^Kh K {r). 



(41) 
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Figure 1. Numerical solution of n x vs x to the Field (1959b) 
problem for a unit source with a Doppler profile. Shown are the 
analytic solution (solid), the FPA solutions to second order (short- 
dashed) and third (long-dashed) order, and the FPDA solutions of 
Rybicki & Dell' Antonio (1994) (dotted short-dashed) and Rybicki 
(2006) (dotted long-dashed). Also shown is a truncated second 
order FPA solution (dotted line; see text). The effect of atomic 
recoil is included for a medium at T = 10 K. Top left panel: The 
even part of n x at the times, from bottom to top, r = 0.1, 1, 10, 
100 and 1000. Only the second and third order FPA solutions are 
shown at r = 100 and 1000 for clarity. In the absence of recoils, 
the even part is identical to the full solution. Top right panel: 
The odd part of n x . The odd part vanishes in the absence of 
atomic recoil. Bottom left panel: The specific temperature T„(x) 
at t — 1000. Bottom right panel: The evolution of the mean 
harmonic temperature (T n )n in terms of the heating efficiency. 
The points correspond to the approximation Eq. (19) applied to 
the Rybicki (2006) solution. 

It is noted that even for an even source in x, cosmological 
expansion will introduce an imaginary part to h K . It follows 
from Eq. (39) that energy exchange will occur even without 
atomic recoil. This results in a cooling term as the photon 
distribution is redshifted through the line-profile (cf, Chen 
& Miralda-Escude 2004). It is also noted that for a linear 
inflow of matter, 7 changes sign and produces a heating term 
as photons are blueshifted through the line-profile. 



4 RESULTS 

Field (1959b) derived an analytic solution to the problem 
of the scattering of resonant photons in a homogeneous 
and isotropic medium with a Doppler profile source term. 
The problem is not particularly well-suited to the Fokker- 
Planck approximation because of the rapidly varying source 
function, but it does provide a test of the approximation 
and its convergence. Solutions to Eq. (35) are found with 
S{x) = exp(-x 2 )/7T 1/2 for four cases: the second and third 
order FPA equations, and the FPDA equations of Rybicki 



& Dell' Antonio (1994) and Rybicki (2006). The results for 
cold hydrogen gas at T = 10K are shown in Figure 1. The 
solutions are separated into their even and odd parts, as, 
for a source even in x, the effect of atomic recoil appears 
only in the odd part of the solution, which is proportional 
to e. (The contribution of recoils to the even part is of or- 
der e 2 .) The original Rybicki & DelP Antonio approximation 
best matches the analytic solution near the line-centre. The 
solution of Rybicki (2006) agrees well, but not quite as well 
as the original Rybicki & Dell' Antonio (1994) form. 

By contrast, the second order Fokker-Planck approxi- 
mation derived here shows appreciable features. The sub- 
traction of a 3 /6 from a/2 in Eq. (38) in the Rybicki (2006) 
formulation plays the critical role of eliminating the oscil- 
latory behaviour. This corresponds to the highest (fourth 
order) derivative of 4>v{x) in Eq. (34). Truncating the equa- 
tion at the third derivative of 4>v(x) by subtracting this term 
from the second order solution eliminates the oscillatory be- 
haviour. Accordingly, a truncated second order FPA (TFPA) 
equation is solved with this term removed. The term, how- 
ever, is removed only from the equation for the real part of 
n K , not the imaginary, for the reason explained below. Al- 
though the (full) third order Fokker-Planck solution agrees 
more closely with the analytic solution than the second or- 
der, the oscillatory behaviour is still present. 

Field (1959b) demonstrated that when atomic recoil is 
included, the radiation field reaches statistical equilibrium 
with the matter near the line-centre, in the sense that for 
e\x\ << 1, the radiation field asymptotically approaches the 
form n x ~ n ( x e) [l-2ex] ~ n x e) exp[-ft(i/-i/)/fcT] for r >> 1. 
It follows from Eq. (21) that the colour temperature at line- 
centre is T„(0) — T. Since the analytic correction to the 
profile due to recoil is proportional to x, the extent over 
which the profile assumes this form is unclear. Field (1959b) 
argued it should apply within a flattened core with the time- 
asymptotic value n x {j) ~ r/[2(log r) 1 ^ 2 ] at line-centre, at 
frequencies smaller than a critical frequency x c — (log r) 1 ^ 2 . 
In Figure 1, the odd part of the analytic solution is plotted as 
n x °^ — n£ e ' [exp(— 2ex) — exp(2ex)]/2, noting that it is strictly 
applicable only for e\x\ « 1. The slow logarithmic growth 
of the critical frequency suggests the radiation field may 
reach statistical equilibrium outside the Doppler core only 
very slowly. Moreover, the value of the colour temperature 
depends also on how flat the core is away from line-centre. 
This suggests the heating rate may not become vanishingly 
small for many scattering times at line-centre (r >> 1). 

Figure 1 shows that the frequency-specific colour tem- 
perature T n {u) approaches the gas temperature at line- 
centre only for the Rybicki & Dell' Antonio (1994) and Ry- 
bicki (2006) approximations. Even for these cases, T n (y) de- 
viates considerably from the gas temperature well within the 
core of the line profile. In the other cases, the colour tem- 
perature takes on mostly very small values, both positive 
and negative. The lesson would appear to be that this tem- 
perature has little physical meaning, except possibly very 
near the line-centre. The colour temperature in the Fokker- 
Planck approximation, however, appears to converge to the 
equilibrium temperature at line-centre only very slowly with 
increasing order. Much away from line-centre, it is not even 
clear the colour temperature should converge to the gas tem- 
perature, except for the extreme blackbody limit. 

As discussed in the previous section, both the heat- 
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ing rate and the Wouthuysen-Field efficiency depend on the 
mean harmonic temperature {T„)h- Weighting by the fine 
profile produces a mean harmonic temperature that evolves 
toward the matter temperature for large r, as shown in Fig- 
ure 1. 

ft is noted that while subtracting a 3 /6 from a/2 in the 
equation for the real part of h K in Eq. (38) removes the 
oscillations in the second order FPA solution, subtracting 
this term from the equation for the imaginary part as well 
produces a solution that results in a rapid, and unphysi- 
cal, decline in (T„)h to values well below the gas temper- 
ature. Subtracting the term only from the real part yields 
{T„)h — > 9.8 K ~ T by r £ 100, and stays constant at this 
value up to r = 1000. The light temperature has reasonably 
converged to the matter temperature within the error of 
the approximation. The truncated version of the second or- 
der FPA equation thus appears a viable approximation for 
including the effects of atomic recoil while recovering the 
correct solution to the radiative transfer equation, in spite 
of the lack of a mathematical justification. 

The form of T„(v) in the different approximations has 
a curious property. While it agrees nearly exactly with the 
gas temperature in the Rybicki & Dell' Antonio (1994) and 
Rybicki (2006) Fokker-Planck diffusion approximations at 
line centre, this is no longer true in the full Fokker-Planck 
approximations. In these cases, the cancellations among the 
negative and positive values of T n (v) ensure the mean har- 
monic average colour temperature {T n )n converges to the 
gas temperature. For the Rybicki & Dell' Antonio (1994), 
Rybicki (2006) and truncated second order Fokker-Planck 
approximations, the convergence is rapid, with the efficiency 
factor scaling roughly like 1 — T/{T„)h ~ r -1 ^ 2 for r >> 1. 
For the other cases, the convergence starts rapidly, but by 
r > 100 levels off at efficiency values of 5 — 10% by r = 1000. 
It is unclear how large r must be before the efficiency de- 
clines to 1% and heating essentially shuts off. Going to third 
order indeed suggests the time is longer than given by the 
second order Fokker-Planck approximation. The integration 
method used is too inefficient to integrate to much longer 
times or to higher orders, so that it is not possible to answer 
this question here. 

The points in the lower right hand panel of Figure 1 cor- 
respond to Eq. (19) applied to the Rybicki (2006) solution. 
At early times this results in efficiency factors exceeding 
unity, which is unphysical. At late times, the approximation 
agrees reasonably well with the Rybicki (2006) predicted ef- 
ficiency. 

The case of a constant continuum source S(x) = 1 is 
shown in Figure 2. The gas temperature is T = 10 K, corre- 
sponding to a = 0.015 for the Voigt profile. The second and 
third order FPA solutions again show substantial structure, 
and the efficiencies reach plateaus for r > 100 at values of 
15 — 20% by r — 1000. Again, heating persists at a non- 
negligible level for many scattering times at line centre. For 
the two FPDA and the TFPA solutions, the heating effi- 



ciency continues to decline, again scaling 



jhly 



-1/2 



for t >> 1. The points in the right hand panel correspond to 
the approximation Eq. (19) for the efficiency applied to the 
Rybicki (2006) solution. The values are again unphysically 
high at early times. 
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Figure 2. Numerical solution of n x vs x for a constant con- 
tinuum source. A Voigt scattering profile is adopted with a gas 
temperature T = 10 K, corresponding to a = 0.015. Shown are 
the no-scattering solution (solid), the FPA solutions to second or- 
der (short-dashed) and third order (long-dashed), and the FPDA 
solutions of Rybicki & Dell' Antonio (1994) (dotted short-dashed) 
and Rybicki (2006) (dotted long-dashed). Also shown is a trun- 
cated second order FPA solution (dotted line; see text). The effect 
of atomic recoil is included for a medium at T = 10 K. Left panel: 
The evolution of n x at the times, from bottom to top, r = 0.1, 
1, 10, 100 and 1000. Right panel: The evolution of the mean har- 
monic temperature (T n )jj in terms of the heating efficiency. The 
points correspond to the approximation Eq. (19) applied to the 
Rybicki (2006) solution. 



5 DISCUSSION & CONCLUSIONS 

5.1 Fokker-Planck approximations 

Previous formulations of the time-dependent Fokker-Planck 
approximation to the radiative transfer equation for the 
scattering of resonant photons were based on a Taylor series 
expansion of only the radiation intensity in the scattering in- 
tegral (Unno 1955; Harrington 1973; Rybicki & Dell'Antonio 
1994; Rybicki 2006). The result of the approximation has the 
form of a diffusion equation. In this paper, an equation based 
on a Taylor series expansion of the product of the scattering 
probability and the intensity is obtained, which is the stan- 
dard Fokker-Planck approximation. This permits an easier 
analysis of the improved accuracy of the approximation to 
higher orders. In this paper, the resulting approximation has 
been derived for Case II scattering using the full Voigt pro- 
file, and including the first-order contributions from atomic 
recoil. The approximation correctly conserves particle num- 
bers with well-defined coefficients, unlike the Fokker-Planck 
diffusion approximation. 

Time-dependent solutions are presented to the FPDA 
and FPA equations for the case of a unit Doppler profile 
source (for Case I) and for a constant continuum source. 
An additional term in the second order FPA equation intro- 
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duces artificial features within the cores of the profiles. These 
diminish in the third order, but are still pronounced. Trun- 
cating the second order FPA equation by removing this term 
eliminates the features. In principle, the features should be 
eliminated by going to a sufficiently high order in the Fokker- 
Planck expansion. Such an investigation, however, requires 
an alternative integration method to that used here. 

5.2 Approach to statistical equilibrium 

An expression for the heating rate of the scattering medium 
through atomic recoil is derived. It is shown that the Fokker- 
Planck approximation correctly recovers the full heating rate 
to the accuracy of the approximation. Both the heating rate 
and the effectiveness of the Wouthuysen-Field mechanism 
depend on the mean harmonic colour temperature {T u )h 
computed from the frequency-specific temperture T u {v), 
weighted by the absorption profile. For photons of frequency 
v in statistical equilibrium with the scattering medium at 
temperature T, it has usually been expected that T u {v) con- 
verges to T for large r, where r is the time expressed in 
units of the scattering time at line-centre. This was shown 
to be the case analytically by Field (1959b) for two particu- 
lar problems involving Case I redistribution, at least near the 
line centre. Numerical solutions are found to Field's prob- 
lem for a pure Doppler profile source. Using the FPDA ap- 
proximations of Rybicki & Dell' Antonio (1994) and Rybicki 
(2006), T u (v) ^Tis found at the line centre for r >> 1, but 
away from the line centre T u (v) deviates appreciably from 
T, and even takes on negative values. For the FPA equa- 
tions, |T„(i/)| is generally small compared with T, and takes 
on both negative and positive values. The third order FPA 
equation predicts the larger value at line centre. Presumably 
higher order Fokker-Planck approximations must yield the 
correct temperature at the line-centre. The investigation of 
this question is in progress. 

The more physically-motivated quantity {T u )h is found 
to converge toward T. The rate of convergence, however, is 
sensitive to the form of the approximation. The most rapid 
convergence is found for the truncated Fokker-Planck ap- 
proximation. Whether or not this is an artifact of the form 
of the approximation is unclear. Convergence is also found 
for the Fokker-Planck diffusion approximations. For the full 
Fokker-Planck approximation developed here, the conver- 
gence of {T u )h to T substantially slows for r > 100. Going 
to third order results in an even slower convergence rate. 

5.3 Angular dependence 

The results presented have all assumed an isotropic radia- 
tion field. In situations in which the gas is heated by one or 
more local sources, the incident radiation field may be highly 
anisotropic. One approach is to take angular moments of the 
radiative transfer equation. Neglecting stimulated emission, 
the upper level population in equilibrium is given by 

n u (u)A u i = ni^- <j> J dv R(i/',ti;v, vl)I v , (r, t, n'),(42) 

and the emissivity becomes 

r\ v = hvni^- (p j dv' R(v , n'; v, n)/„/(r, t, ft'). (43) 



The angular coupling between the redistribution function 
and the intensity precludes a simple form for the energy 
transfer rate. A solution to the radiative transfer equation 
along these lines is sketched by Rybicki & DelP Antonio 
(1994). 

5.4 Astrophysical implications 

The Wouthuysen-Field mechanism for coupling the spin 
temperature to the colour temperature of resonant photons 
has been discussed in a variety of contexts, including the 
IGM, intergalactic gas clouds like Damped Lyman-Alpha 
Absorbers, gas on the peripheries of galaxies, and in clouds 
near quasars, with applications not only to H, but D and 
3 He+ as well (eg, Field 1958; Deguchi & Watson 1985). The 
mechanism may be particularly complicated for 3 He + be- 
cause of the Lya radiation from 4 He + and the coincidence 
of the 3 He + Lya transition with a Bowen fluorescence line of 
O in . The reader is referred to Deguchi & Watson (1985) for 
a discussion of these and other complicating effects. The dis- 
cussion in this section is confined to more general remarks. 

The scattering times at line-centre for 1 H, 2 D and 3 He + , 
normalised to the hydrogen number density riH (in cm~ 3 ) 
using the Big Bang Nucleosynthesis constrained abundances 
of 2 D and 3 He from Buries, Nollett & Turner (2001), are 

t sca = 3.2n^T 1/2 s ^H, 
= 7.3 x 10 4 n^T 1/2 s ; 2 H, 

= LOxlO^gV^s ; 3 He + . (44) 

These timescales are typically short for neutral gas. For ion- 
ized gas, however, as for instance in Lyman Limit Systems, 
the rates may become quite long, particularly for D and 
3 He + . For instance, in a cloud with nm ~ 10~ 5 cm~ 3 at 
T = 10 4 K, the timescale is about 10 12 s for these isotopes. 
If it takes 10-100 scattering times for the Wouthuysen-Field 
mechanism to bring the colour temperature (T u )h to the 
gas temperature, the equilibrium timescale approaches that 
of the lifetime of massive stars. If such objects, or their 
H ii regions, are the sources of the Lya photons, it may be 
that while the colour temperature at the H i resonance has 
come into equilibrium, that of D i has not, in which case 
an erroneous D/H ratio may be inferred from the compari- 
son of radio measurements of the two hyperfine transitions. 
Conversely, if the cosmic D/H ratio is accepted, such ob- 
servations may provide an estimate of the time since the 
radiation sources turned on. Even for hydrogen, the colour 
temperature may not reach equilibrium if the source varies 
on a timescale that does not greatly exceed the scattering 
time. 

Cosmological expansion will generally play a negligible 
role in the approach to statistical equilibrium of the radia- 
tion field. This is because the smallness of the Sobolev pa- 
rameter. Rybicki & Dell'Antonio (1994) show that, within 
the damping wings, expansion negligibly affects the photon 
distribution until a time r « (a/7 4 ) 1 / 3 >> 1 for typical 
values of a and 7 expected for the IGM. 

Expansion, however, may play a role in more local situ- 
ations, such as the outflow from a star or an active galactic 
nucleus (AGN). One possible effect is a freeze-out of {T u )h- 
For the sake of demonstrating the principle, a computation 
with a = 0.015 and 7 = 0.1 is show in Figure 3 for cold hy- 
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Figure 3. Numerical solution of n x vs x for a constant continuum 
source in an expanding medium with 7 = 0.1 A Voigt scattering 
profile is adopted with a gas temperature T = 10 K, correspond- 
ing to a = 0.015. The results for the truncated FPA solution 
are shown in the left panel at the times r = 0.1 (dotted), 1 (dot 
short-dashed), 10 (dot long-dashed), 100 (short-dashed) and 1000 
(long-dashed). The photons are increasingly redshiftcd with time. 
The light temperature (T n )u is shown in the right panel for the 
truncated second order FPA equation (solid line) and the Rybicki 
(2006) FPDA equation (dashed line). 



Figure 4. Numerical solution of n x vs x for a constant continuum 
source in an expanding medium with 7 = 0.01. A Voigt scattering 
profile is adopted with a gas temperature T = 2 X 10 4 K, corre- 
sponding to a = 3.3 X 10~ 4 . The results for the truncated FPA 
solution arc shown in the left panel at the times r = 0.1 (dotted), 
1 (dot short-dashed), 10 (dot long-dashed), 100 (short-dashed) 
and 1000 (long-dashed). The photons are increasingly redshiftcd 
with time. The light temperature (T n )jj is shown in the right 
panel for the truncated second order FPA equation (solid line) 
and the Rybicki (2006) FPDA equation (dashed line). 



drogen gas at T = 10K. The photon distribution is skewed 
to the red by the expansion. Both the truncated second or- 
der FPA solution and the Rybicki (2006) FPDA solution 
predict a freeze-out of {T u )h above the gas temperature for 
r > 100, but at different values. The radiation density pro- 
files differ slightly near line centre, with the dip at x — in 
Figure 3 not as deep in the FPDA solution. 

For hot ionised hydrogen gas, the expansion may re- 
sult in a light temperature well below the matter temper- 
ature. This is illustrated in Figure 4 for a — 3.3 x 10~ 4 
and 7 = 0.01 for gas at T = 2 x 10 4 K. The photon dis- 
tribution is strongly skewed toward the red by r = 1000, 
with a sharp cut-off in the blue similar to the Chugai (1980) 
steady-state limiting solution for the radiative transfer in 
the profile wings of photons from a monochromatic source 
at line-centre in an expanding medium. Both the truncated 
second order FPA solution and the Rybicki (2006) FPDA 
solution predict (T U ) H < T for r > 60. The Rybicki (2006) 
FPDA predicts negative values of {T u )h for r £ 50. It should 
be noted, however, that at the assumed temperature the re- 
quirement 2e << b/c is only marginally met. The full Ry- 
bicki (2006) formulation includes an extra term of order b/c 
which is not incorporated into the formalism here. This term 
may eliminate the negative values, although it is noted that 
negative values are still found in a similar integration using 
T = 1000K. 

The possible heating role played by atomic recoil was 



originally proposed by MMR as a possible means of heating 
the IGM at high redshifts before it is re-ionised. Given the 
rapidity with which {T u )h —> T, this now appears unlikely 
(cf, Chen & Miralda-Escude 2004; Rybicki 2006), particu- 
larly if the predictions using the Rybicki & Dell' Antonio 
(1994), Rybicki (2006) or truncated second order Fokker- 
Planck approximations are most accurate. The prediction 
of the FPA equations, however, is a levelling off of the ef- 
ficiency to values of 10-20% for large r. This suggests that 
heating may persist for r >> 1. How long is still unclear. 
The delay may also be an artifact of the artificial oscillatory 
behaviour of the solutions. 

The heating rate of MMR at full efficiency may be ex- 
pressed as 

Gmmr = hU " hv mPi. (45) 
m a c 2 

The amount of energy per atom in the lower level that may 
deposited before the efficiency becomes vanishingly small 
may be expressed as 

AE = f(G M MR/ni)t sca = f h »° h (46) 

m a c z m 

where n VQ = f dv(f>v{v)n v . For the hydrogen 
Lye* transition, hvo/m a c 2 ~ 10 -8 and hvo Si 10 5 K. Conse- 
quently a very large ratio of photons within a Doppler width 
of the resonance to neutral hydrogen atoms is typically re- 
quired for an appreciable increase in temperature. Expressed 
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in terms of the average photon occupation number N(fo) at 
line-centre, 

AS = /- 5-- TT^hvo 



m a c 2 c n ; Ag 



8 x 10" 



A/~M 
io- 21 



K, 



(47) 



where the expression has been evaluated for hydrogen 
Lyo? in the last line, and normalised to JV(u ) = 10 -21 , 
corresponding to a mean intensity of J„ = 2.2 x 
1CP 22 ergcm~ 2 s _1 Hz -1 sr _1 . Thus in a normal intergalac- 
tic environment, the effect is completely negligible even prior 
to the reionization epoch. 

There may be some circumstances, however, where the 
effect may contribute non-negligibly to the total heating. For 
instance, for blackbody radiation at T = 10 4 K, M{yo) ~ 
exp[— hva/kT] « 8 x 10~ 6 . For a somewhat dilute black- 
body not in thermal equilibrium with the matter, recoils 
may provide a significant amount of heating. Whether or 
not it would be competitive with other processes depends 
on the particular situation. 

At the opposite extreme, cold tenuous gas near a lumi- 
nous source such as a young star or active galactic nucleus, if 
shielded from ionizing radiation, could be heated by several 
degrees or more. 

Another situation in which recoil heating may play a 
role is in an H n region. According to the on-the-spot ap- 
proximation, the photon occupation number at line centre 
due to diffuse recombination radiation is (Osterbrock 1974), 

1,2 \ 3 / 2 i 

— 7lHI,(48) 



2nm e kT 



ram 



4xw~ 22 t; 



Xhi 



where T 4 is T in units of 10 4 K, Xhi is the neutral hydrogen 
fraction, and m e is the mass of an electron. Using this in 
Eq. (47) gives 



AE » 3 x 10" 



(£) 



t: 



Xhi 



K. 



(49) 



It is useful to express this in terms of a source of ionising 
photons S for a Lya optical depth thi through a region of 
size R. Then 

AE « SSsoTHr'.Rpc -1 ^ 1 K, (50) 

where 6*50 is in units of 10 50 s" 1 and R pc is in units of 
parsecs. The scattering of diffuse Lya photons may pro- 
vide a significant source of heat in the development of an 
H 11 region near an ultraluminous star or an AGN. 
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APPENDIX A: MOMENTS OF THE 
REDISTRIBUTION FUNCTION 

The moments {Q)<p(v) of the redistribution probability 
W(v,Q) are computed here. These results generalize those 
of Rybicki & Dell'Antonio (1994) to explicitly take into ac- 
count atomic recoil. Accordingly, the results are derived here 
from first principles. 

Consider the absorption of a photon of frequency v trav- 
elling in direction n and re-emission of a photon of frequency 
v' travelling in direction fi' by an atom travelling with veloc- 
ity v in the laboratory frame. It is useful to project v onto 
an orthogonal basis given by the average and difference of 
the photons directions (Hummer 1962): 



ni = 7+(n + n), n 2 =7-(ft- 



n 3 = n x n',(Al) 



where j± = 1/[2(1 ± /it)] 1 / 2 and fi — n ■ n'. The correspond- 
ing components of the atom velocity are (vi , vi , V3) . To low- 
est order in v/c and huo/rn a c 2 , where m a is the mass of the 
recoiling atom, the frequency shift on scattering is given by 

hi/ 



(ft' 



r(l-M) 



(A2) 



Allowing for radiation damping in the upper level, the 
absorption profile is given by the Lorentz profile in the rest- 
frame of the atom 



m = 



S/tt 



[{v - ^o) 2 + <5 2 ] ' 



(A3) 



where vq is the line-centre frequency and <5 = r/47r, where 
r is the radiative damping width of the upper state. In ex- 
pression Eq. (A3), the implicit approximation v/uo ~ 1 is 
made, except for the frequency shift in the denominator, and 
so does not recover the correct behaviour far from line cen- 
tre, such as the Rayleigh limit (Peebles 1969; Peebles 1993). 
Since this paper is concerned solely with energy transfer 
resulting from resonant scattering, this is an adequate ap- 
proximation. Radiation heat exchange at other frequencies 
must be dealt with separately. 
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It is assumed that the atomic velocites are distributed 
as a Maxwellian with temperature T according to 

P(u)d 3 u = 7r~ 3/2 exp(-u 2 )d 3 «, (A4) 

where u = v/6, and b = (2kT/m a ) 1 '' 2 is the Doppler pa- 
rameter. It is useful to introduce the Doppler width ui = 
uob/c, frequency offset x — [y — vo)/u, recoil parameter 
e = (hfo/rn a c 2 )c/b = hvo/(2kTm a c 2 ) 1 ^' 2 , and normalised 
damping constant a = 5/u). The frequency shift may then 
be expressed as x — x' — 7Z 1 U2 + e(l — /i). Then the probabil- 
ity that an incident photon of frequency shift x is scattered 
to x' = x — q is given by 



w(x,q,n) = 



dui / dv,2 exp(— u\ — u 2 ) (A5) 

; J — oo 

£[g ~ 71^2 - e(l - Ap] 

a 2 + [x - 7+tti(l — 7-U 2 (l - At)] 2 ' 

neglecting terms of order xb/c in the denominator. 

It is noted that the form of w(x, q, n) satisfies the asym- 
metry relation Eq. (8) in the limit hvo » kT and v jv ~ 1 
(to order xb/c). The reverse scattering of x — > x — x — q is 
x' = x — q — » x' — (—<?). The corresponding relation between 
w(x,q,Li) and w(x — q,—q,fi). Expressing w(x — q,—q,n) 
according to Eq. (A5) and making the change of variable 
U2 = U2 — 2"f-q shows after straightforward manipulations 
that 



w(x-q,-q,n) = e 2tq w(x,q, li) 



— exp 



kT 



w(x,q,fi). (A6) 



The angle- averaged coefficients {Q n )ip(y) may be ex- 
pressed in terms of w(x, q, fx) according to 



(Q n M 

where 
a n (x,n) 



\ 1 71-1 f 



dfia n (x,Li), 



dqq n w(x,q,Li). 



(A7) 



(A8) 



Following Rybicki & Dell' Antonio (1994), such integrals are 
most easily evaluated in Fourier space. From Eq. (A5), the 
Fourier transform of w(x, q, li) is 



du2 exp(— u 2 — u\) 



w(k,X,h) = J dx J dqexp(ixn) exp(iqA) 
a r . f 

^L dui J- 

S[g - r )Z 1 u 2 - e(l - fi)] 



a 2 + [x - 7+ui(l + n) — 7_U2(1 - n)] 
= exp 



1 2 

— a\n\ — -k +i\e(l — fi) 



--A K (l-/i)--A 2 (l- M ) . (A9) 
The Fourier transform of a„(x,Li) may be expressed as 



noting that factors of q n correspond to n derivatives with 
respect to iX in Fourier space. The first five values of o„(k, fx) 
are then 



a (K,Li) = 

ai(K,Li) — 

a 2 (K,Li) = 

& 3 (k,li) = 



exp 



I 1 1 a' 

-OK -K 

4 

1 



e(! - A*) + - At) 



x exp 



1 1 2 

— OK K 

4 



(1 - A*) + (1 - M) 2 (« 



1 2\ 

4 K ) 



1 1 2 

-OK -K 

4 



x exp 



3(1 -A«) (e+2 iK 



(All) 



(A12) 



(A13) 



i(l - A*) 3 (6ek 2 +iK 3 ) 



x exp 



11 1 2 

-OK K 

4 



(A14) 



0,4 (K, //) 



3(1- M ) 2 +6(1- M ) 3 ( 



1 2 

^eK -K 

4 



+ ^(1-A0 4 (-8ieK 3 +K 4 ) 



16 

x exp 



I 1 2 

—OK K 

1 4 



(A15) 



The n — term corresponds to the full absorption pro- 
file <p(v). Since for Case II redistribution, this corresponds to 
the Voigt profile, the notation (pv(f) is adopted for ip{v). Its 
dimensionless form will be denoted as <f>v{x). The Fourier 
transform of 4>v(x) is given by oo(k, li) in Equation (All). 
The first five moments of Q are then 



(Q°) = l; 

(Q) = £ q,-l^ ^^M ; 

(Q 2 > = ^-l^ rflog^W 
3 Of 

d 2 ip v {v) 



1 4 1 
-UJ 



3 ¥>v(f) d^ 2 



(A16) 
(A17) 

(A18) 



(Q 3 ) = 4 e a, 3 -2a, 4 dl0g ,f v(t/) +^ 5 \ /^ u) 
dv 2 </'v(i / ) dv i 



i^e 1 

4 ¥'v(t / ) d^ 3 



(Q 4 ) = 4a; 4 - 126c; 

8 7 1 



5 dlog^v(!/) 



dv 

' 3 <pv(v) 



(A19) 

d 2 tp v {v) 



pv(v) dv 2 



5 ¥'v(i / ) d^ 3 



1 



d 4 v?y(^) 



5 H^viy) dv A 



(A20) 



9™ 

On(«,A*) = Q(jxyl W(K,A,M) 



(A10) 
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